function P_SAAD=Pclogit_SAAD_decsunk(Y_SAAD,X_SAAD,Z_SAAD,b_SAAD)

choice=2:13;
ns=numel(choice);

Y=sortrows(Y_SAAD,[1,2]);
Y=Y(:,3:end);
Y=Y-1;  %!!!### UPDATED on 06/03/2020: Y MUST be from 1 to J!!!

X=sortrows(X_SAAD,[1,2]);
X=[ones(size(X,1),1),X(:,3:end)];

Z_SAAD=sortrows(Z_SAAD,[1,2,3]);
Z_SAAD(Z_SAAD(:,3)==1,:)=[];
Z=zeros(size(Y,1),size(Z_SAAD,2),ns);
for i=4:size(Z_SAAD,2)
    Z(:,i,:)=reshape(Z_SAAD(:,i),ns,[])';
end
Z(:,1:3,:)=[];

baseAlt=1;

P_SAAD = pclogit(b_SAAD,Y,X,Z,baseAlt);


